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, ' Abstract 

^ ' We present a general framework to design Godunov-type schemes for multidimensional ideal mag- 

, netohydrodynamic (MHD) systems, having the divergence-free relation and the related properties of 

the magnetic field B as built-in conditions. Our approach mostly relies on the Constrained Transport 
' - ^ ' (CT) discretization technique for the magnetic field components, originally developed for the linear 

^-^ , induction equation, which assures [V • B]„„m — and its preservation in time to within machine 

' accuracy in a finite-volume setting. We show that the CT formalism, when fully exploited, can be 

«j I used as a general guideline to design the reconstruction procedures of the B vector field, to adapt 

■ standard upwind procedures for the momentum and energy equations, avoiding the onset of numerical 

monopoles of 0(1) size, and to formulate approximate Riemann solvers for the induction equation. 
This general framework will be named here Upwind Constrained Transport (UCT). To demonstrate 
^ ■ the versatility of our method, we apply it to a variety of schemes, which are finally validated numeri- 

5— ( ' cally and compared: a novel implementation for the MHD case of the second order Roe-type positive 

, scheme by Liu and Lax (J. Comp. Fluid Dynam. 5, 133, 1996), and both the second and third order 

versions of a central-type MHD scheme presented by Londrillo and Del Zanna (Astrophys. J. 530, 
^ I 508, 2000), where the basic UCT strategies have been first outlined. 

X' 

H ■ 1 Introduction 

In extending Godunov-type conservative schemes designed for Euler equations of gas-dynamics to the 
system of (ideal) magnetohydrodynamics (MHD), in the multidimensional case a main problem arises 
on how to represent the solenoidal structure of the magnetic field vector B and on how to formulate 
reconstruction procedures and (approximate) Riemann solvers sharing consistency with this property. In 
the last years a number of works have focused on this specific problem and many different approaches have 
been proposed. A wide class of (second order) numerical schemes for regular grids have been analyzed and 
compared by Toth PP, while contributions covering also higher order schemes, adaptive mesh refinements 
(AMR) and unstructured grids are in rapid development. 

Since we are mainly interested here to analyze methodological aspects, we propose a broad classification 
of the published contributions on this specific topic into two main groups: 

1. Schemes based on standard upwind procedures (henceforth SUP) designed for Euler equations, where 
also magnetic field components are discretized at cell centers as the other fluid variables. Since in 
this case the approximated [V • B]„„„i based on central derivatives may have a non- vanishing size, 
different strategies to control or prevent the accumulation in time of related spurious numerical 
effects (usually referred to as numerical monopoles) have been proposed. 

• A first method, suggested by is to add an elliptic (Poisson) equation to recover the 
solenoidal property at each time-step. In reference PP this procedure has been named pro- 
jection scheme and is currently widely adopted (see [3] for a high order WENO scheme). 
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• In the scheme introduced by Powell 0] (see also [S]), the numerical [V • B]„um quantity is not 
forced to vanish; the MHD system is reformulated by adding new source terms proportional to 
this variable in order to recover the original MHD system in non-conservative form. Moreover, 
the classical seven-mode Riemann wave fans have been enlarged to eight modes. In this 
modified system, upwinding is applied to all magnetic field components and hence also to the 
component i3„ across a discontinuity surface. 

• In a more recent work |S| , in order to preserve both the conservative form and the hyperbolic 
structure of the MHD system, a new time dependent wave equation is introduced to damp 
and/or to transport away the non-zero [V • B]„„m contributions. 

2. In the second group we include schemes which take advantage of the so-called Constrained Transport 
(CT) method by Evans and Hawley [7] (originally suggested for the evolution of the induction 
equation in the linear approximation). It is a main feature of this method to introduce staggered 
discretizations of magnetic and electric vector fields in the induction equation. In fact, by using 
these staggered values to approximate the relevant first derivatives, [V • B]„„m = in the initial 
conditions and its exact preservation in time result. The problem here is on how to apply this 
formalism in a Godunov-type scheme for the full MHD system. 

• Most of the published works combine the above CT discretization with the SUP cell centered 
discretization by introducing different empirical recipes (e.g. [S|) EIj QHIi E])- However, these 
procedures result in a sort of hybrid schemes and the problem of numerical monopoles is still 
left open, in our opinion. 

• In our previous work 12 (LD from now on) we have proposed numerical procedures to take 
advantage of the specific CT discretization benefits, and hence the [V • B]„„m = condition, 
even in the reconstruction steps and in the approximate Riemann solvers. The same method 
has been then applied to relativistic MHD 

The goal of the present paper is twofold. First, by adding analytical arguments to the approach 
outlined in LD, we propose a method to construct and then to characterize a class of numerical schemes. 
Second, we present implementations of a variety of different schemes, to demonstrate the versatility and 
self-consistency of the method. 

Regarding the first goal, our main concern is here to select a set of properties, some of them common 
to the Euler system and other specific of MHD equations, which in our opinion should be preserved in 
the numerical discretization. In this way, it is then possible to envisage Godunov-type schemes for MHD 
having: (a) the divergence-free condition as an exact built-in property, (b) reconstruction and upwind 
procedures consistent with this property. Since the CT formalism comes out to be the necessary starting 
point to achieve this result, our framework will be named here Upwind Constrained Transport (UCT) 
method. 

As a novel numerical application we then propose the UCT implementation of the positive scheme 
by Liu and Lax (|21j ISDi a second order Roe- type scheme which proves to be accurate and robust. 
Numerical validation will be finally presented for several standard two-dimensional test problems, where 
the results of the new MHD positive scheme are compared with central-type schemes as proposed in LD, 
extended here to more accurate central-upwind two-speed approximate Riemann solvers, and tested in 
its second and third order implementations. 

This paper is organized as follows. In the next sub-section we propose and discuss some general 
conditions as guidelines for numerical modeling. The main ingredients to formulate general UCT-based 
Godunov-type schemes for MHD systems, i.e. the discretization form, the proper reconstruction proce- 
dures and the approximate Riemann solvers, are presented in Sect. 2. In Sect. 3 we specify the method 
to the positive and central MHD schemes, which will be finally tested and compared in Sect. 4. 

1.1 Conservation laws and consistency demands for numerical MHD 

The MHD system has a peculiar form and cannot be simply reduced to a set of conservation laws for 
scalar variables, as the Euler equations. In fact, if the specific structure of spatial differential operators is 
taken into account, it is more properly represented by the set of the following two coupled sub-systems: 

^+V.f(w)=0, (1) 
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— + VxE(w) = 0, (2) 
equipped with the non-evolutionary constraint on the B vector field 

V • B = 0, (3) 

which, once satisfied for initial conditions, is analytically preserved in time by Eq. 1^)1. 

The set of equations evolves in time the five-component array of scalar functions u = [u'(x, t)]"^, 
I = 1,2, ... ,5, while the set (O evolves the vector field B — [Bi{x,t)]'^, i = x,y,z. The overall set of 
dependent variables are henceforth represented by the eight-component array w — [u, B]-^. The first 
array contains the conservative fluid variables u = [p,qi,e]'^, where p is the mass density, qi = pvi are 
the momentum components, Vi are the fluid velocity components, and e = p/(7 ~ 1) + pv^ /2 -\- /2 
is the total energy density for a perfect gas equation of state, where p is the kinetic pressure and 7 is 
the adiabatic index. The corresponding flux vector components = [flY" , I = 1,2, ... ,5 are given by 
fj = \<li, Mij, Hi]'^ , i, j = X, y, z, with the momentum flux tensor defined by Mij = ViQj + liSij — BiBj 
and the energy fiux components defined by Hi = Ui(e -I- 11) — _Bi(v • B), where 11 = p + B^ /2 is the 
total pressure. In sub-system which is the induction equation for the magnetic field vector B, the 
corresponding flux is simply given by the electric field vector E = — v x B, where the assumption of a 
perfect conducting plasma (ideal MHD) has been implicitly assumed. 

As for the Euler equations, the system (|1I2|I has to be supplied with entropy functions S = S'(w) 
satisfying the condition 

dS 

— + V-Fs(w)<0, (4) 
ot 

which allows to identify, among discontinuous solutions of the MHD system, the (physically) admissible 
ones. The existence of entropy functions (in fact S — ~ps, where s cx log{pp^'^) is the physical entropy 
per unit mass) is also related to the hyperbolic structure of the MHD equations. For smooth solutions, 
the system (|1I2|I can be put in the non-conservative (quasi-linear) form 

^ + [J(w) . V]w = 0, (5) 

where J = (Ji), i = x,y,z, and each Jj is the Jacobian matrix of the eight-component flux array 
[fl,Ei]'^ with respect to the w variables. It is a well known property that any linear combination 
J(w, k) = fci Ji(w), for real ki numbers, and then also each Ji matrix, is hyperbolic at any reference 
state w. Moreover, as for the Euler equations (see ^^), the (positive) Hessian matrix S'w,w acts as 
similarity transform to make all Ji symmetrizable. 

To underline differences and analogies of the MHD system with respect to the reference Euler system 
which may have relevance for numerical modeling, some remarks are in order: 

• The u array contains scalar variables and the corresponding flux derivatives are expressed by the 
div = [V-] conservative operator. Sub-system has then the same formal structure of the Euler 
system for gas-dynamics. At surface elements where discontinuities take place, this conservation 
form leads to the usual Rankine-Hugoniot relations. On the other hand, the B(x,t) vector is anti- 
symmetric (an axial vector), components Bi are pseudo-scalars, and the corresponding evolution 
operator is given by the anti-symmetric curl = [V x •] derivative. The conservative form is now 
expressed by the scalar condition Q (magnetic flux conservation) and by the V x E flux derivatives 
(conservation along a closed contour). Discontinuous solutions satisfy jump relations just for the 
tangential components Bt = B x n, where n indicates the normal direction, whereas the normal 
field component i?„ = B • n is continuous. The Rankine-Hugoniot relations, once supplied with 
an appropriate entropy law, allow to identify the physically correct discontinuous solutions. It is 
apparent that magnetic discontinuities and the related entropy constraint do not involve the parallel 
Bn component. 

• It follows that smoothness properties of MHD variables are also different. Scalar components w'(x, t) 
may develop discontinuous solutions along any space direction and can be then represented on the 
space of piecewise continuous functions. The vector field B(x, t) has more elaborate properties, 
since the divergence- free condition entails the B(x) field maps piecewise differentiable (and then 
continuous) field lines. The conservation law given by Eq. |(5J is then essential to preserve in time 
condition Q and to assure the smoothness properties of the magnetic field. 
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• The divergence-free condition enters implicitly in the MHD momentum and energy conservative 
equations. This can also be expressed by realizing that the Maxwell tensor T — IB^ /2 — BB in 
the momentum flux has to satisfy 

B • (V • T) = 0, (6) 
in order to recover the correct Lorentz force in non-conservative form. 

• Finally, the divergence- free condition allows to represent the B(x, t) field via a vector potential 
A(x, t), defined by B = V x A and by the gauge condition V • A = 0, which assures the uniqueness 
of this representation. The new evolution equation is now 

f+E = 0. ,7) 

The above relations and the induction equation ||2Jl, together with the condition E • B = valid for 
ideal MHD, imply an added conservation law for the magnetic helicity H — / (A • B)da;'', carrying 
informations on the topology of magnetic field lines. 

When looking at (finite-dimensional) numerical approximations, a main problem is that no rigorous 
results on convergence are available, even for the Euler system. In this case, however, by taking advantage 
of theoretical achievements, like the Lax-Wendroff theorem |17| . heuristic guidelines are usually adopted 
in order to: 

• retain the conservative form of the original equations in the discretized system; 

• assure consistency, in the sense that the approximations of the flux functions and of the differential 
operators have to recover the exact ones as the spatial and temporal grid sizes go to zero; 

• assure non-oscillatory (or even monotonicity preserving) numerical representation of discontinuous 
data; 

• assure consistency with the entropy law, in a way the numerical viscosity induced by the upwind 
differentiation is compatible with Eq. |0J (see |18p: 

• assure stability of the numerical solution. 

As already anticipated in the Introduction, the main issue addressed here is to select a set of additional 
requirements for the MHD system which should assure that the specific properties of the magnetic field 
enter as built-in conditions of a numerical scheme. We propose the following: 

• the discretized first derivatives diBi entering the V • B definition are consistent approximations; 

• for initial divergence-free fields the approximated derivatives satisfy [V • B]„„„i — exactly; 

• divergence-free initial conditions are preserved exactly in time by the discretized induction equation. 

We then suggest the following definition: a numerical scheme is consistent with the specific properties 
of the MHD system if all above conditions are fulfilled. This definition, together with the guidelines for 
Euler equations, will enable us to identify and construct a class of Godunov-type schemes for MHD, later 
referred to as UCT-based schemes. 

In this framework, as for the Euler equations, a finite volume setting provides a sufficiently general 
starting point. Here we concentrate only on algorithms for regular structured grids, even if the generality 
of the method allows to extend some basic procedures also to adaptive mesh refinements (AMR, |^) and 
to unstructured grids. In particular, De Sterck 20 has developed a general CT formalism for unstructured 
triangular grids, named MUCT, where rigorous geometrical arguments have been considered to support 
this approach. 
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2 The UCT method to design Godunov-type schemes for MHD 



2.1 Discretization step: finite- volume formalism 

In a finite-volume setting, the 3-D computational domain Q is first subdivided in Cartesian cells C, with 
volume V, side sizes hi, i = x,y,z and faces given by the oriented surface elements S^, i — x,y,z, 
where ± denotes the sign of face normals. For each face Si, we then denote as {L^,L^), j,k ^ i its 
oriented sides. At this level of the analysis no indexing on a grid is needed, thus allowing to extend the 
formalism also to a non-uniform partition of the domain, as required in grid refinement techniques. 
In the following, a semi-discrete finite-volume approach will be employed, thus only space averages will 
be considered and the time dependency will be left for further integration, for example via standard 
Runge-Kutta algorithms |21| . 

A conservative discretization of sub-system is accomplished, as usual for Euler system, by integrat- 
ing each scalar equation on the volume element V of each cell C. By application of the Gauss theorem, 
one has then 

^"(^) + E^(f^-fn = 0, (8) 

i 

where 

m ^^J^ "(^' f f = ^ /^^ f« [w(x, t)]dS (9) 

denote respectively volume averages of each scalar component over the cell C, and f ^ are flux values av- 
eraged on cell faces . We note that fi fluxes are represented as exact point values in the (non averaged) 
parallel coordinate i, and the corresponding differences in provide the averaged flux derivatives. 

In the case of sub-system two different approaches can be pursued. In the SUP approach, magnetic 
fleld components are discretized by volume averages Bi as other scalar variables u', and the electric fleld 
components by face averages Ek, as fluxes. We have then 

^B.W+^e.,,,fe^(E+-Ej=0, (10) 

where Sij^k is the Levi-Civita symbol and ± here refers to faces normal to the j direction. 

On the other side, in the CT formalism a discretization preserving the original (vector anti-symmetry) 
property is accomplished by a surface integration on a cell face followed by the application of the Stokes 
theorem on the line contour of that face. This leads to 

^6.(i)+5]e.,,fc^(E+-Ej = 0, (11) 

where now 

Ut) = ^ / B,{^,t)dS, Et = TL: f Ek[M^,t)]dL (12) 

are respectively the staggered discretized magnetic field variables, defined as integrals over the cell face 
Si (we retain the formalism of non-capital bi components to indicate staggered values to conform with 

other authors) , while Ef. are now line-averaged electric field components along face edges , where the 
orientation depends on the normal to the face under consideration (see Fig. Here the magnetic field 
components are represented as (normalized) magnetic fluxes, thus bi are exact point values in the parallel 
coordinate i (as fi) while Ek are point values with respect to their orthogonal coordinates {i,j) ^ k. 

By using this staggered discretization, which yields a couple of normalized fluxes 6^ defined at Sf faces 
for each direction i, it is now possible to represent the volume average of the (parallel) first derivatives 
diBi{x). Therefore, for divergence-free initial conditions Q, we have 

E/J-(^^^-^n = 0, (13) 

i 

which will be preserved in time algebraically by the induction equation (|ll|l . 

At this general level, anything is exact. Approximations (in space) arise when all MHD variables, 
starting from the discretized values, have to be reconstructed at the cell faces where fluxes are defined as 
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point values. However, even at this preliminary step, differences in tlie magnetic field representations of 
Eq. Hll|) with respect to Eq. (|f 0|) have relevance: 

• When primary data for the magnetic field are the staggered bf components, on each cell C one has 
at disposal two independent sets of data. A first consequence is that no reconstruction is needed 
to evaluate (at a second order approximation) these variables as argument of the corresponding 
fi fluxes at a Si face. Moreover, these staggered data carry informations both on the volume 
averaged (or centered) values Bi and on the first derivative along the parallel coordinate. In fact, 
(foj — 6j )/hi provides a second order (and then consistent) approximation of the diBi first derivative 
in point-wise (or finite difference) sense inside each cell. 

Related to the above is the property that each component provides a continuous sampling across 
the corresponding Si cell face. This follows from the definition in Eq. (|12|) and by taking into 
account that a divergence- free B(x) field entails a continuous elemental flux across a discontinuity 
surface (see |12| for details). In this way, the continuity property of the _B„ normal component in 
the Rankine-Hugoniot jump relations retains a consistent representation in a finite volume setting. 

Staggered components can also be defined, in a fully equivalent way, by using the vector potential 
A. In fact, by face-averaging the defining condition B = V x A, one has 

6. = ^e..,,fc^(A+-:4J, (14) 

still assuring the divergence-free condition in the form H13|l . If needed, a time evolution for the 
numerical (line averaged) Ak components, consistent with Eq. Q, can also be formulated: 

^Ak{t)+Ek = 0, (15) 
at 

preserving in time the representation of Eq. (|I4(I . Here, the numerical flux Ek are precisely the 
same as in Eq. 

• On the other hand, when primary data are represented by Bi volume averages, as in Eq. Hl()|l . only 
one numerical value per cell is available. The averaging procedure, which is now applied also along 
the i coordinate, while still assuring that Bi are consistent approximations of Bi point values, it also 
entails a loss of direct information on the point values at the Si cell faces and on the correspondent 
parallel derivatives. Moreover, interpolation procedures are not sufficient, for discontinuous data, to 
recover these informations (related to the divergence-free property) and then some added argument, 
procedure, or constraint would be necessary, in our opinion, to avoid inconsistent approximations. 
In the following section we shall provide a more detailed analysis on this problem. 

We can then conclude here that a finite volume CT formulation of the induction equation satisfy 
the the general consistency demands presented in Sect. 1.1, provided the diBi are approximated using 
staggered bi data. 

In spite a CT based discretization for MHD is a longstanding well known framework for induction 
equation, its application to the whole MHD system is yet a matter of debate. In fact, it is yet a persistent 
viewpoint in numerical community that staggered collocation of magnetic field may be useful only to 
express the divergence-free relation, being otherwise not well suited for upwind formulation in Godunov- 
type schemes. In contrast with this viewpoint, we propose here to construct and test numerical schemes 
explicitly based on the CT discretization and on the related properties detailed above. 

2.2 Reconstruction step: scalar vs divergence-free fields 

In the following, we specialize to a Cartesian partition of the computational domain Q. For grid indexing 
1 < J < Nx, 1 < fc < Ny, I < m < Nz, a generic cell Cj,k,m is defined as 

Cj.k,m = [2^^-1/2,2:^ + 1/2] X [yk-l/2,yk+l/2] X [z,n-l/2,Z„i+i/2], (16) 

where each fractional index labels a cell interface (say here with < j < Nx), while a cell 

center has coordinates {xj,yk,Zm), with Xj = {xj^i/2 +Xj^i/2)/2, yk = {Vk-1/2 + Vk+i/2)/'2, and Zm = 
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+ •2^m+i/2)/2- For simplicity, now we assume a uniform partition, along all directions, so that 
{h^ = a;j-+i/2 - Xj_i/2, hy = yk+1/2 - yk-1/2, = 2:„+i/2 - 2^-1/2) are the constant sizes for all 
cells. Under these settings, the primary volume-averaged array of fluid variables will be indicated as 
Uj,fe,m- Surface-averaged staggered magnetic field components, defined at cell interfaces, will be indicated 
as {hx)j+i/2,k,rm (by) j,k+i /2.7m ^ud (^z ) j,fc,m+i/2 • The Same notation holds then for face-centered flux 
components f^, i = x,y, z, while the edge-centered Ek fields have indexing (£'z)j+i/2,fe+i/2,m: and so forth 
for other components. Finally, the divided differences for functions located at cell interfaces introduced 
in the previous section, e.g. in Eqs. 0, lfTT|l . and will be denoted here as Di, i = x,y,z. For a 

generic 1-D scalar function / located at inter-cell points Xj±i/2, is then given by 

[D,{f)], - ^[A./],; [Ax/], = /,+i/2 - f,-i/2- (17) 

In higher order Godunov-type schemes, a one-dimensional scalar variable u{x), represented by cell- 
centered data {uj}, is first reconstructed as approximated point values u{x) inside any cell Cj and up to 
the interior cell faces Xj±i/2, where fluxes fx{u) have to be evaluated. This is accomplished by using local 
polynomials u{x): a) consistent with the cell averages values Uj] b) having monotone or non-oscillatory 
properties. In a second order approximation, one has the linear fit 

Uj{x) — Uj + Dx{u){x ~ Xj), (18) 

where the non-oscillatory derivative Dx(u) is usually constructed using slope limiters. In the simplest 
case it is defined as ^ 

[Dx{u)]j = —mm{[Axu]j+i/2, [Axu]j-i/2), (19) 

where [A2;u];-|-i/2 = ~ui (I — j, j — 1) and where mm{a, b) denotes the usual two-point MinMod (MM) 
algorithm. More elaborate limiters can likewise be constructed, using the Van Leer j22| monotonicity 
constraint or TVD properties. For higher order schemes, ENO-based procedures have been developed 
(see for a review) assuring a non-oscillatory reconstruction under weaker monotonicity constraints. 

The reconstruction in Eq. 1)1811 extends up to the interior face points Xj±i/2 providing a left approx- 
imation on the Sx face and a right approximation on the S~ face, along the indicated coordinate, as 
needed for flux computation. At a given cell interface, say dX x — Xjj^i/2-, in cases where a jump of 
size Aa;M = 0(1) occurs, the estimated slope coefficients \Dx\j and using minmod limiters both 

vanish and the reconstructed uix) variable is represented by piecewise constant (first order) interpolants 
Uj — Uj and Uj+i = Wj+i, respectively. 

For a multidimensional function u(x) (in the uniform grid defined above), a tensor-product represen- 
tation with one-dimensional interpolants on each coordinate is usually adopted. For unstructured grids 
more elaborate procedures are required, but the basic ingredients (consistency with the cell averaged 
data and non-oscillatory constraints) still hold. In regular grids, the resulting second order approxima- 
tion u{x) « u{x) on each cell Cj_k.m takes then the form 

t2(x) = u + Dx{u){x - Xj) + Dy(u){y - yk) + Dz{u){z - z,n), (20) 

where all quantities are implicitly calculated at cell center and each Di (u) , is the non-oscillatory 1-D first 
derivative defined in Eq. H19|) . 

When dealing with the magnetic field B(x), a different approach is needed, in general, to take into ac- 
count the vector structure and the specific smoothness properties already quoted in the previous sections. 
Starting with face-averaged data bi, for i = x,y, z, satisfying the divergence-free condition 1)13(1 

Dxibx) + Dyiby) + D,{h) = 0, (21) 

the problem of representing each Bi (x) field along the proper parallel coordinate inside a cell is already 
solved at the linear level, since the slope Di(bi) is at disposal. Instead, reconstruction is needed along 
the face (orthogonal) coordinates, where the field is sampled by averaged ^x{x)]k,m data. 
For the -Ba;(x) function one has then 

S,(x) =Bx + Dx(bx){x - Xj) + Dy{Bx){y ~ yk) + Dz{Bx){z - z,„). (22) 

Here, the cell-centered Bx values result by noticing that Bx{x, yk, Zm) is the the unique linear interpolant 
of a continuous function with data at the Xj±i/2 points. We have then 
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{Bx)j,k,m = -j[{bx)j + l/2 + (&x)j-l/2]/c,m, (23) 

which provides an approximation of the cell averages and of the Bx point values: 

Bx^Bx + 0{hl)^Bx + 0{hl). (24) 

In Eq. (|22|l it is evident how staggered and cell-centered fields work differently. Along the parallel x 
coordinate, the Bx{x) function is entirely defined by the {hx)j±i/2,k,m data, providing the field values at 
the cell interfaces, the approximated first derivative [Dx{bx)\j.k,im and the approximated cell-averaged 
value {Bx)j,k,m- Reconstruction is needed instead to evaluate the Bx left-right values at the orthogonal 
cell faces [S^, ^t): where Bx may have discontinuities and where, correspondingly, bx and cell-centered 
values Bx behave as the other u scalar variables. 

A second property related to the divergence-free conditions is that the Bx{x, •) function in Eq. (|22|l 
maps a continuous function with first derivative Dx{bx) which may be discontinuous. To evaluate the 
jump size, one considers the second difference A^b^, = Ax^xbx centered at the proper interface point 
{xj^i/2, yk, Zk)- By taking into account the divergence-free relation 1)21(1 and the commutativity property 
of the two-point difference operators, one has 

-^Ap^ = --^AyAxby ~ -^A.Axh, (25) 

x y z 

where at least one of the differences Axby or Axbz has size 0(1), by definition. 

Following the same procedure the (By,Bz) components can be represented, on the same cell, by the 
relations 

Byi^) = By + Dx{By){x ^ X ,) + Dy(^y){y - yk) + Dz{By){z - Zra), (26) 

B,(x) - + bx{B,){x -x,) + Dy{B,){y - yk) + D,(b,){z - z„), (27) 

and remarks made above on different behaviors of staggered (hy,bz) and cell-centered {By,Bz) values, 
respectively, and on the different smoothness properties depending on the involved coordinates, apply. 

We remark here that this apparent duality in the reconstruction procedure appears to be fully consis- 
tent with the physical duality of the MHD Rankine-Hugoniot relations. In fact, using a local characteristic 
decomposition of the (u, B) variables, for a specified direction, say x, only the (u. By, Bz) variables par- 
ticipate to the Riemann wave fan and may contribute, then, to the discontinuous characteristic modes. 
On the other hand, the continuous bx variable has no role in the related upwinding procedures. 

In upwind schemes, the reconstruction step presented above has relevance not only to recover face 
centered values, but also to approximate the B(x) field at any point inside a cell, as it is required in 
schemes adopting grid refinements (AMR) or multi-grid procedures. In this context, divergence-free 
interpolants similar to those derived here, even if based on quite different arguments, have been proposed 
in PH. A related work has presented a detailed analysis to show that, under linear interpolation 
based on staggered data, conservative properties and the divergence-free relation can be preserved in 
cell-subcells refinement procedures. 

To summarize, it follows from this analysis that a numerical divergence-free magnetic field can be 
represented in an unambiguous way by using bi staggered values (or equivalently the related numerical 
vector potential) as primary data. In this representation, second order approximated first derivatives are 
consistent and non-oscillatory (no cell crossing is needed) and the variable [V • B]„„m = '^i[Di{bi)\j,k,m 
in Eq. (|21|l results to be exactly zero inside any point of Cj^k.m- At the same time, the reconstruction 
procedures at cell interfaces, where variables are discontinuous, provide definite rules on how to formulate 
upwind differentiation. 

2.2.1 Reconstruction and central derivatives in non-CT schemes 

A CT-based formalism helps also to analyze the reconstruction problem in the SUP framework, where 
only cell averaged Bi values are at disposal. By restricting to the Bx{x, •) function, a linear interpolant 
reads 

Bxix, •) = [Bx]j + [Cx(Bx)]j{x - xj) + .. (28) 

where Cx{Bx) denotes a (by now unspecified) consistent approximation of the first derivative. For a 
piecewise differentiable function, at a point Xj^i/2 this interpolant provides in general two values, as 
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left and right approximations, like for all other variables. By imposing there the additional continuity 
condition, specific to magnetic field components, one has 

[Bx]j + 2f^x[Cx{Bx)]j = [Bx]j+i — -hx[Cx{Bx)]j+i, (29) 

and since Taylor expansion for is not applicable across a discontinuous interface, this assures only an 
implicit way to express the numerical derivative in terms of the cell centered data. This suggests that some 
added condition, like the non-oscillatory constraint, should be required even for a;-wise interpolations, 
but how to recover the divergence-free condition remains here an open question. 

We are not aware of any SUP-based scheme where this problem has been properly addressed. It is a 
common practice, instead, to use for flux computations the mid-point average 

{ix)j + l/2,k,m = ■^[{Bx)j.k,7n + {Bx)j+l.k,7n], (30) 

resulting in a 0{hx) approximation. In turn, this entails an approximation for the first derivative given 
by the central difference 

[B>i''\Bx)]j,k,7n = -^[{Bx)j+l,k,m - {Bx)j-l,k,7n]- (31) 

However, when the relevant (a;j_i, Xj+i) stencil includes a discontinuity, say at Xj_^_i/2, by using the 
relation H29|l with Cx — Dx(hx) we have 

[Di'\Bx)],,k.m = [Dx{hx)]j,k.,n + 0(1), (32) 

the 0(1) term resulting from the first derivative jump as estimated in H25|l . Using the same argument to 
the other {By,Bz) components, a final [V • B]„„m = 0(1) results. This is well documented in numerical 
experience where a central difference is used for discontinuous functions (the Gibbs pathology) . 



2.2.2 Extension to higher order 

The first derivative discontinuity of a divergence-free field has also relevance to extend a CT-based recon- 
struction to higher orders (r > 3) of spatial accuracy. In Godunov-type schemes, higher order reconstruc- 
tions for scalar variables are usually provided by ENO-based interpolants, like Weighted-ENO (WENO, 
US' and references therein) and Convex-ENO (CENO, f57j). For regular grids, these interpolants, once 
defined for one-dimensional variables, can be extended to higher dimensions by a tensor-product repre- 
sentation. 

For the magnetic field more elaborate procedures are needed, however, by taking into account that 
non-oscillatory derivatives along different directions are non-commutative. A general strategy we propose 
is to reconstruct first the vector potential components Ai in the usual way, by taking advantage that these 
are scalar variables, and then to define the bi point values using the basic B = V x A relation. However, 
using the vector potential alone is not sufficient to guarantee a divergence-free relation. A crucial step 
to make this procedure effective is to approximate the V x A derivatives by consistent two-point, fixed- 
stencil, high order finite differences. 

As an example, in the Appendix we report the third order implementation of the CENO reconstruction 
procedures. 

2.3 Upwind step: Roe-type approximate Riemann solvers 

Using the grid notation of the Sect. 2.2, the MHD Eqs. JSJ and (|1I|I take on the form 

^[u(t)] + ^A(fO = 0, (33) 

i 

and 

^[S^(i)]+E^M,fc^.(^^)=0, (34) 
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the first set being centered at cell nodes and the second at cell interfaces. The overall system has now 
to be evaluated using some approximate Riemann solver, the same for all flux functions and for all 
Cartesian components, at a time. Here we make reference to Roe-type schemes, allowing a full resolution 
of the characteristic MHD modes, whereas the so-called central and central-upwind schemes, which avoid 
spectral decomposition, will be briefly treated in Sect. 3.2 

Let us then first specialize to the 1-D flux differentiation along the x coordinate and denote with 
Fa; — [F,^]'^ , for s — 1,2, ... ,7, the array of all flux components, defined at the x = Xj^i/2 point for generic 

{y, z) coordinates of the — face. These components are — fx \ for Z = 1, 2, . . . , 5, whereas the 
components entering the induction equation are fJ^-* — —E^ and fJ^-* — Ey. Correspondingly, we denote 
with w^; — [p,qi,e. By, the array of variables which need reconstruction as point values at the Sj, 
face and as [w^{y, z), {y, z)] (East- West, see Fig.|21) the corresponding left-right states. The magnetic 
field component bx{x, •) in the parallel direction i satisfies — by continuity, thus the arguments 
of flux functions Fj;(w) have then to be specified as {w'^,bx), for a = E,W. For short, we denote as 
FS = F,(w^,6,). 

Approximated Riemann solvers based on local linearization technique (see e.g. |2H1: for the MHD 
case) rely on the Roe matrix A^, defined as numerical Jacobian by 

Ff - Ff = i,(w) . (wf - wf ) ^ Ax ■ 5xWx, (35) 

at any {y, z) point of the indicated Sx cell face. As usual, the Ax matrix is evaluated at an appropriate 
intermediate state w = [wa;(w^, wf'), 6^;], and is required to be consistent with the Jacobian matrix 
Jx presented in Sect. 1.1. We notice the Ax matrix has rank seven, as in the pure 1-D case where 
Bx = const holds. In fact, in the multidimensional case the continuity condition 5xbx — plays a similar 
role, implying Bx is locally constant and does not participate to the characteristic wave fan. 

To express the local flux variations in Eq. H35(l in terms of characteristic modes, let us consider the 
spectral decomposition Ax = [RKR~^]x, where, if \s,s = 1,2,..., 7, are the Roe matrix real eigen- 
values, A = diag{As}, and columns of the R matrix are the corresponding right eigenvectors. In this 
representation, let us then split Ax — [Ax\^ -f [A^;]", where the first term contains components with 
Xs > and the second with As < 0, respectively. In this form, one has also \Ax\ — [Ax]'^ — [^k] , where 
\Ax\ — [R\A\R~^]x and |A| = diag{|As|}. By using standard procedures, from the two-point flux function 
Fx{'w^ ,bx) a one-valued, continuous and monotone flux can be selected as upwind state by 

F^ = Ff + [Ax]- ■ Sxv^x = - [Ax]+ ■ 6xWx, (36) 

so that the usual Roe flux formula comes out as 

F^(wf,wr,6,)=F;-*., (37) 

where 

F: ^ i(Fr + Ff ); = • (wf - wf ), (38) 

the first term expressing the smooth component leading to a centered two-point formula in flux differen- 
tiation and the second the Roe-type component coming from the upwind procedure. At a discontinuity 
interface, the latter provides numerical dissipation which needs to be consistent with entropy conditions. 
To that purpose, for shock solutions where As ~ 0, a small amount of added dissipation As — > As +?7s has 
to be introduced as entropy-fix to avoid unphysical behaviors. 

With a straightforward extension, for coordinates i — y, z, the Ai Roe matrices and the upwind fiuxes 
FY{wi,bi) — F* — can likewise be constructed, each function being evaluated at points of the Si 
proper orthogonal face. 

By taking advantage of the Roe formalism, which is based on independent 1-D Ai matrices, where 
characteristic modes are represented locally as planar waves along each direction, it is possible to evaluate, 
at a time, numerical fluxes collocated at different points. In fact, in Eqs. (|33|) the final five-component 

_ (1 — 5) 

numerical fluxes ~ Fj for i = x,y, z are obtained as an average over the proper face Si. This average 
involves only interior face points and then only a characteristic wave fan, the one represented by the Ai 
matrix. On the other hand, the remaining f'^ flux components appearing in the induction equation 
()34|l are defined as point values at the intersections of cell faces, where different characteristic wave fans 
overlap. These flux components can be likewise evaluated by a linear combination of 1-D upwind fluxes 
along the intersecting direction. It is a main feature of the UCT method that this combination follows a 
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proper upwind selection rule, since a same flux component at the same collocation point results to have 
two independent representations in terms of characteristic wave fans. 

As a prototype, we consider the flux, which is defined at the (a;j-|_i/2, yfc+i/2j 2) points where 
faces Sx and Sy intersect. Let then denote as Wj, = [p,qi,e, B^, B^f" the set of variables having a 
representation in terms of the Ay matrix eigenmodes and by Fy = (w^ , w^, the corresponding 
flux array, where now (w^, ) (North-South) denote the left-right states along the y coordinate. At the 

indicated intersection points, = — i^i^-* — Fy^'^ comes out to be a four-state function Ez{yv°"''') where 
a = N, S, b = E, W, since the w argument contains both and Wj, variables (see again Fig.|21). 

The sixth flux component Fx^^ = —E^, defined in Eq. (|37|l and specialized at the y — yk+1/2 point, is 
represented by two independent contributions coming from states (w^, w"^): 

[^^.^(w^)]^^^ - ff " 4>^x, [F^{-^')r = ff - (39) 

where /* — F*^^^ and (px = ^^x^ ■ On the other hand, the Fy^^ — E^ flux function, defined for generic x 
values, is represented by two independent contributions for states (w^,w^) at the x = Xj^if2 point: 

[^^,^(w^)](^^ = ff - [F;/{^'^)r = ff - (40) 

where now /* = Ff^ and 4>y = By taking into account that F^ = —F^'^ at the same 

(cCj_|_i/2, yk+1/2) point, a one- valued numerical flux function having the continuity and upwind properties 
along each direction can be then constructed at the linear level by 



where 



E''A^)=El-<t>y + 4>.. (41) 

K^\[E^^ + El^ + E^^ + Er]- </>, = i(^^^ + ^^^), = + (42) 

and E"^'^ = Ez{^°''^), with a = N, S, b = E,W. The continuity property of this relation implies, 
in particular, that for any possible orientation, namely along a cell face diagonal or along a cell face 
side, the corresponding one-dimensional planar mode is taken into account in proper way, to within 
0(|(5j;Wj;| |(5j^Wj,|), by the upwind combination above. The Ez flux, given in Eq. H41|) for generic z coor- 
dinate of the [Sx, Sy) common side, is finally evaluated by line averaging as {Ez)j+i/2.k+i/2,m numerical 
flux. In a similar way, the other x and y components can then be constructed by a proper combination 
of upwind fluxes along the corresponding orthogonal coordinates and line averaging along the parallel 
coordinate. 

This completes the presentation of the UCT formalism. For later reference, we quote here the flux 
formulas given above, specialized to a first order approximation, which constitutes the building block of 
any Godunov-type scheme, needed for higher order extensions. A piecewise constant reconstruction for 
Wx variables at the x — a;j+i/2 point, needed in Eq. l|S7|l . gives = {wx)j, = {■Wx)j+i, at any 
{y,z) point. Face averaging reduces to the one-point evaluation Fx{'W'x,bx) = Fxi^Wxjbx) and the first 
order numerical flux in Eq. (|37|l , for the five-component f flux array, reads then (over-bars are henceforth 
omitted for brevity): 

i^x)j+l/2,k,m = (f^)i+l/2,fc,m - i^x)f_^^if\,k,m^ (43) 

where, at the indicated {yk,Zm) points, 

(C)j + l/2 = i[f.((w,), + i, {bx), + l/2) + UMj, {bx)j + l/2)], (44) 

and ^ 

[^x]j + l/2 = ^\Ax{Wj + i/2)\ ■ [{'^x)j + l - {'^x)j], (45) 

and so forth for other directions. To the same first order approximation, the Ez flux in Eq. 141|) can 
be calculated at (a;j+i/2, J/fe+1/27 ^m) edge points (the Zm centering will be assumed implicitly in the 
following). By specializing the Ez = f* ~ —f* ~ —{vxBy — VyBx) argument variables, the smooth term 
results as: 

{E*z)j+l/2,k+l/2 = -7;[iVxby)j + l + {Vxby)j]k+l/2 + 7;[{Vybx)k+l + {Vybx)k]j+l/2, (46) 
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where {vx)i,k+i/2 = ■^[ivx)k + {vx)k+i]i, for I = j,j + 1, and {vy)j+i/2,i = ^[{vy)j + {vx)j+i]h for 
I = k,k + l. It is worth noticing that in Eq. (|46|l . the E* term contains the {bx, by) staggered components 
and the resulting four-states combination at the (a:^j+i/2i yfe+1/2) point cannot be reduced simply to an 
interpolation or averaging form based on the four cell centered values of the argument. On the other 
hand, in the dissipative Roe-type fluxes, centering at y = yk+1/2 of the (j)x term comes out as a two-point 
average in the orthogonal coordinate and correspondingly for the (py term, so that 

{4>x)j+l/2,k+l/2 = l^[{'t'x)k + {4'x)k+l\] + l/2, {4'y)j + l/2,k+l/2 = ■^^{<Pv)j + )i+l] fe+1/2 • (47) 

On the computational side, this form is also economical and of easy implementation, since the four 
contributes there involved can be derived from the ^x and $y fluxes already worked out for the fluid 
variables. 

2.4 On the problem of numerical monopoles 

We consider now some main differences of the present approach with other schemes, by focusing on 
the problem of numerical monopoles. These unwanted effects may arise when the ^iDi(ti) term for 
momentum equations in Eq. H33|l fails to recover the proper [J x B]„tjm Lorentz force in the original non- 
conservative form, or, equivalently, when the orthogonality condition ® is not satisfied with sufficiently 
high accuracy. 

This problem has been analyzed in details by j3U| . in a general discretization setting and with no 
particular reference to specific upwind differentiations. Here we follow some of his arguments and no- 
tations to show the behavior of various classes of MHD schemes in comparison with our UCT method. 
We specialize to second order fiux differentiation, as it is implemented in most of the schemes in the 
literature. 

By restricting for sake of simplicity to a two-dimensional configuration, the Maxwell stresses Tij , i,j = 
x,y, entering the momentum flux components are given by 

Tx.x = -Ty^y = -{By - BI), Tx,y = Ty^x = -BxBy, (48) 

evaluated at the proper interface points, that is ixj±i/2,yk) for T^^x and Ty^^^ (x^ , ?/j.-|-x/2) for Ty^y and 
Tx,y We notice that these terms refer specifically to the /* smooth part of the relevant fiux components, 
since Roe-matrix contributions give only a diffusive term of the Lorentz force. For fiux differences related 
to the Qx momentum, one has then the algebraic relations 

DxiTx.x) = tix{By)Dx{By) ~ fix{Bx)Dx{Bx), Dy{Ty^x) = -fiy{By)Dy{Bx) - ily{Bx)Dy[By), (49) 

all being centered on a {xj,yk) point, where Dr^ and Dy are the usual divided differences of Eq. 117|) . 
and where for components a = Bx,By we define the averages fix{o-) = io-j+i/2.k + ^Jj-i/2,fc)/2 and 
^^yi^) — {0'i,k+i/2 + ctj,fe-i/2)/2. These two-point averaging on the cell center give Hxio-) — o-j,k + 0{h^) 
and Hy{a) = aj^k + 0{hy), for second order numerical fluxes. By summing the two differences defined 
above, one has (let h = hx = hy) 

DxiTx.x) + Dy{Ty^x) = -Bx[Dx{Bx) + Dy{By)] + By[Dx{By) - Dy{Bx)] + 0{h''), (50) 

where the second term on the right hand side provides the numerical approximation of the x component 
of the Lorentz force. By considering then the correspondent qy momentum flux, the approximation to 
the orthogonality condition finally becomes: 

[B • (V • T)]„„™ = -{Bl + Bl)[V ■ B]„„™ + 0{e), (51) 

which is satisfied to within the truncation error if the numerical divergence-free relation holds (at least) 
to the same accuracy order. For smooth fiows this is clearly true in any discretization of first derivatives. 
But when discontinuities are present, this term can be even of 0(1) size. 
In particular: 

• In schemes where (-Bx)j+i/2.fc and (Sj^)j.fc+i/2 are reconstructed as two-point average of the cor- 
responding cell-centered data, as discussed in Sect. 2.2.1, the resulting [V • B]„„m in Eq. H51|l is 
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expressed by the central differences oi (Bi), giving contributions of 0{h^) size in smooth regions 
and of 0(1) size near discontinuous interfaces. 

In Powell's eight-wave scheme, the source terms introduced in the momentum equation given by 
B([V • B]„„m) yield a subtraction of these monopoles, which is another way to express the Lorentz 
force in its non-conservative form. 

In projection schemes, a [V • B]„„,„ = condition, yet based on cell centered data, is enforced at 
each time step as an added constraint, so that monopoles, which arise when flux derivatives are 
computed, are prevented to grow in time. 

In previous CT-based schemes (e.g. [HI, [3], CODi where staggered variables are actually at disposal, 
only cell-centered averaged data Sj, derived as in Eq. I|23|) . are then employed in flux computations. 
In that case the resulting [V • 'Bi]num, term in Eq. (|51|l is still expressed via derivatives based on 
central differencing as in Eq. H31|l . and numerical monopoles are now produced. These can be 
evaluated now exactly by noticing that 

D^^\b,) ^ D^{b,) + hlR,, R,^^^Alb, (52) 

and by similar expressions for the other {By,Bz) components. In this relation the residual h^Rx 
is 0(1) at points where the first derivative is discontinuous, as can be seen from Eq. H25|l . The 
final estimate gives [V • B]„„m = J2i hfRi = since the residuals Ri, i = x,y, z, being related 
to the transverse jumps, do not cancel out, in general. 

• The class of UCT schemes proposed here prevents the onset of monopoles, since bi staggered data 
enter directly in the corresponding /* flux components, and [V • b]„„,„ — '^.i[Di{bi)\ — results in 
Eq. ISH). 

We finally notice that, in order to avoid numerical monopoles, flux derivatives have to be computed 
at the same time-stepping level, since time-splitting techniques prevent exact cancellation of [V • 'S\num 
terms. 



3 Examples of application of the UCT method 

In the present section, the numerical strategies outlined in Sect. 2 will be applied to a couple of existing 
Godunov-type schemes, originally designed for fluid dynamics. We have chosen two completely different 
schemes, a classical Roe-type scheme based on field-by-field limiting along characteristics, and a simple 
central-type scheme adopting a two-speed upwind fiux with component- wise limiting in the reconstruction 
algorithms. Both schemes are proposed here in the semi-discrete form, appropriate for our UCT method, 
and then integrated using TVD Runge-Kutta time discretizations of the appropriate order j21| . 

3.1 Roe-type: the positive scheme 

We present here the MHD implementation of a second order fiux-limited scheme proposed by Liu and 
Lax ('|14|. |15p. This approach allows for an easy formulation for multidimensional hyperbolic systems 
and can then represent a well behaved alternative to standard TVD-based schemes (see 221) • Moreover, 
these authors have introduced a new positivity principle, which is more appropriate for multidimensional 
systems, to which TVD do not apply. This stability principle relies, in particular, on the symmetrizable 
form of the system under investigation and is then well suited also for MHD equations. 

Flux-limited schemes are constructed as a proper combination of an accurate second order smooth 
numerical flux, for example the centered approximation F* = (Fj -|-Fj+i)/2 or a Lax- Wendroff term, and 
a first order diffusive flux of the form F''*** = F* — in a way that when the flow is smooth F*^ — F* 
and when discontinuities are present F*^ = F'*'"*. In a symbolic way, this combination can be expressed 
as F^ = F* — (/ — L)#'^\ where L is a diagonal operator whose entries are flux limiter functions (psi^s) 
acting, in general, on characteristic modes and assuring (j)si(^s) = 1 for smooth modes and 4>s(Qs) = 
otherwise. 

In particular, in this positive scheme by Liu and Lax, two first order dissipative schemes are combined 
with different flux limiters: a first (least dissipative) flux of Roe-type and a second (more dissipative) 
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flux of LLF-type, the latter acting also as entropy-fix for the former. In the MHD implementation, the 
Roe matrices Ai(w) are constructed in the space of primitive variables, with a simple two-point average 
to define the state w and with the eigenvector formulation given in 31 to remove degeneracies. 
The resulting formula, for each cartesian fiux component, reads then: 

F7+1/2 = ^[FK) + FK+i)] - ^i?[cROE |A|(/ -L) + CLLF «(/ - M)]R-\w,+, - w,), (53) 

where the two coefficients combining the two dissipative terms are chosen such as croe + cllf = 1 , Af is 
a diagonal term made up of smooth minmod-type limiters that satisfy < (j>s{ds) < 1, < 4>s{0s)/ds < 1, 
L is a diagonal term whose sharper limiters satisfy < 4>s{ds) < 2, < (j>s{Ss)/Ss < 2, for a resulting 
maximum CFL number of 0.5 (for more technical aspects, the reader is referred to the cited works). 

In our experience, when applied to the MHD system by following the UCT strategies described 
in Sect. 2, this scheme results to be a robust and accurate flux-limited scheme, with almost no extra 
computational effort, if compared to standard multidimensional TVD schemes. 

3.2 Central-type schemes 

In the fluid-dynamics community, schemes adopting simple one or two-speed numerical fluxes with 
component-wise reconstruction (no characteristic decomposition is thus required), are often referred to 
as central or central-upwind schemes, since it has been proved these approximate Riemann solvers come 
out in Roc-type schemes by some form of central averaging over the Riemann characteristic wave fan (see 
|32j . and 33 for the latest developments). For most applications, even in shock dominated flows, these 
schemes give yet satisfactory results and are surely more economical than Roe-type schemes, even when 
third or higher order reconstruction is applied. 

In the present paper, two UCT implementations of central-upwind schemes based on the Harten-Lax- 
Van Leer HLL |34j two-speed flux and component-wise reconstruction are considered: 

1. the simplest second order HLL-UCT implementation is performed by reconstructing all (potentially 
discontinuous) variables by a Monotonized Centered (MC) Van Leer 22 limiter. Flux upwinding is 
then achieved by a HLL flux formula, as presented below for the MHD system. Time integration is 
finally given by the second order Runge-Kutta scheme. This second order central-type scheme will 
be here referred to as MC-HLL-UCT; 

2. A third order central scheme, based on the local Lax-Friedrichs (LLF) flux and Liu and Osher 
Convex-ENO (CENG) scheme has been already presented and tested in our first paper LD. 
The choice of CENG has been favored by the following interesting features: first, the high order 
reconstruction algorithm is able to reduce itself to minmod-type limiters (MM or MC) at discon- 
tinuities, thus strongly reducing spurious oscillations typical of component-wise reconstruction via 
ENG interpolants; then, a formulation based on point values, rather than cell averages, allows the 
use of purely one-dimensional reconstruction routines. The same algorithms are used here in a 
central-upwind scheme now equipped with an HLL fiux formula, and this new third order central- 
type scheme will be named CENG-HLL-UCT. For the interested reader, we report in the Appendix 
the main computational steps used. 

In the following, we define the HLL two-speed formulas appropriate for MHD. The HLL upwind flux 
for the fluid components, say fc at {xj^i/2,yk, Zm), may be written as 

^ «+ff + a-ff -«+a-(u^-u^) ^ ^^^^ 
at + ax 

where as usual = f2;(w", bx), a ~ E, W, and 

a± =max{0,±A±(wf ,6,),±A±(wf,5.)}, (55) 

and similarly for the other y and z components. Here, to avoid the definition of an intermediate Roe-type 
state, we have chosen to calculate the dissipative a^ terms by taking the maximum eigenvalues (in MHD 
systems related to the fast magneto-sonic speeds: = Vx i: cl) between left and right states. 

Notice that, by rearranging the terms, it is still possible to rewrite Eq. H54() in the form of Sect. 2.3 
f* — 4>, so that the discussion on [V • B]„„m of Sect. 2.4 applies unchanged. Gn the other hand, we can 
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not use the same composition rules outlined in Sect. 2.3 to derive Roe-type magnetic fluxes. The correct 
form of the upwind flux comes out now by averaging over the two overlapping x and y Riemann wave 
fans at the (x^+i/a, y/c+1/2) edge, thus 

-ta+E^^ + a+a-E?^ + aZa+E^^ + aZa-E?^ 



E 



U _ ^^x ^^y -"^z > "a; "i/ -^z > "a; 



{at + ax ){ay + ay ) 



- -$^ib'x - + - bf), (56) 

ay + ay al + ax 

where the a^ and a^ at edge center points should be calculated by taking the maximum characteristic 
speeds (in absolute value) among the four reconstructed states, whereas for sake of efficiency we actually 
consider the maximum over the two neighboring inter-cell points, where these speeds are already at 
disposal from fluid fluxes computation. 

In spite of a different form, Eq. 1)56(1 share with Eq. ((41|l all the required upwind properties, that is 
to be a four-state formula and to reduce to the correct 1-D flux for shocks aligned with x, y, or x — y 
diagonals. The usual settings are instead recovered for LLF fluxes, which are the particular cases of the 
corresponding HLL ones for a^ ^ a^ = ax and a+ = cny ~ cty, where now each single a^ = max{a^, a~} 
speed defines a symmetrical (central) averaged Riemann fan. 

It is interesting to notice that the upwind fluxes of Eq. 1(54(1 generalize to MHD the formulas of the [HS! 
central scheme for the Euler system. Moreover, Eq. ((56() coincides with that defined in the same paper 
for Hamilton- Jacobi scalar equations, which is to be expected since each component of the induction 
equation, for a given velocity field and expressed in terms of the magnetic potential, has exactly the form 
of such equations. 



4 Numerical results on test problems 

We present now a standard set of numerical examples to assess accuracy, stability and effective divergence- 
free properties of our UCT schemes presented above. In the following tables and plots we shall refer to 
these three schemes as POSITIVE-UCT (the flux-hmited Roe-type scheme of Sect. 3.1), MC-HLL-UCT, 
and CENO-HLL-UCT (respectively the second and third order central-upwind schemes described in 
Sect. 3.2). Most of the choices in the following tests have been inspired by those in T, although a direct 
comparison of UCT with the other methods for MHD proposed or reported is difficult to achieve, due 
to the use of different underlying schemes, and it will be not our main concern here. Comparisons will 
be made instead among our three UCT schemes, and also with respect to their corresponding non-UCT 
counterparts, i.e. when magnetic field components are discretized and evolved exactly as the other fluid 
variables (the SUP framework), with the same underlying scheme. We shall refer to these Euler- like 
versions as base schemes (BS). This will enable us to appreciate the effective improvements introduced 
by the use of our UCT method. 

In the following sub-sections we shall provide quantitative measures of the divergence-free properties 
in our schemes. To distinguish between the two types of numerical representations discussed in Sect. 2.2.1 
and 2.4, here we shall deflne, for each cell Cj,k,m, the two quantities 

[V • B]„„„. = [V • b]„„„ = ^^(^»)' (57) 

i i 

where o'f^ are the central derivatives defined in Eq. 1(31(1 . applied to the cell-centered derived data Bi, 
and Di are the usual two-point divided differences deflned in Eq. 1(17(1 . As already discussed, while 
[V • b]„tim will be zero to within machine accuracy for second order UCT schemes, and arbitrarily small 
(see the Appendix) for CENO-HLL-UCT (it is simply not deflned for the BS counterparts), the [V • B]„t(„j 
variable may have 0(1) jumps at discontinuities. Recall that in UCT schemes the onset of monopoles is 
actually measured by the [V • b]„„„i variable, rather than by [V • B]„„„i like in other CT and non-CT 
methods for MHD where cell-centered flelds are employed in fluxes. Notice that, in spite of the use of 
the staggered magnetic fleld b in the initial conditions and in the computations, the output data will be 
referred to the cell-centered interpolated fleld B. 

To obtain more accurate results in the various test problems, schemes may be tuned by choosing each 
time the most appropriate slope limiter in the reconstruction routine. In all simulations, a CFL number 
of 0.5 and an adiabatic coefflcient of 7 = 5/3 are used, unless differently specified. 
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4.1 A convergence test: the oblique CP Alfven wave 

A propagating circularly polarized (CP) Alfven wave is a well known exact non-linear solution of the 
multidimensional MHD system, and therefore it is often used to measure accuracy of a numerical ap- 
proximation. We notice that this test involves only smooth solutions, thus problems related to the 
divergence- free condition are not expected to arise here. Following jl], we consider a CP wave prop- 
agating on the {x, y) Cartesian plane at an angle a — 30° relative of the x axis. Periodic boundary 
conditions can be applied if the computational box is defined by < a; < 1/ cos a, and < y < 1/sina. 
Let the coordinate along the propagation direction be ^ = a; cos a -\- y cos a and i] = y cos a — x sin a be 
the coordinate along the transverse direction, then the initial values of wave variables may be given as 
Vri = Brj = Asin(27r^) and Vz = Bz = Acos(27r^), where A measures the wave amplitude. The constant 
parallel components are defined as = and = 1, along with a uniform density p = 1 and pressure 
p = 0.1. The chosen values correspond to a wave period T = 1, to a propagation Alfvenic speed = 1, 
and to a sound speed Cs = \/tpIp — 0.4. 

To check numerical accuracy, the evolved solution w(a;, y, t) at a time t,nax = nT, for a given number 
n of periods, is usually compared to the initial condition w{x,y,0) and the difference Sw evaluated in 
some norm is then tabulated at different grid resolutions. To compare with the Toth results, we also 
adopt the Li norm to evaluate the relative errors and the averaged d values are measured only by taking 
into account the transverse vector components w^, Vz, Sjj, and Bz- However, since we have noticed that 
convergence was not precisely achieved with the values A = 0.1 and n = 5 suggested by Toth, especially 
for our most accurate third order central scheme, we have decided to use here the safer values A = 0.01 
and n — 1. We think that the reason is due to compressible effects: in spite of the fact that a CP Alfven 
wave is an exact solution, regardless of its amplitude A, this kind of wave is known to be subject to the so 
called parametric decay instability (see |35| and references therein). This is due to non- linear wave- wave 
interactions that, via coupling to compressive modes, lead to wave distortion and decay, so that when 
this happens we are no longer comparing with the true solution. 

The results are reported in Table for all our three UCT schemes, with resolutions ranging from 8^ 
to 128^. In our highest accurate scheme CENO-HLL we use the smoother MM limiter, while for the 
second order schemes the sharper MC limiter is employed, otherwise results at low resolution are not 
quite satisfactory, and a small value of cllf = 0.01 is used in POSITIVE. In cases where only smooth 
fields are involved, UCT schemes may perform worst than their BS counterparts, due to the additional 
interpolations needed to recover the cell-centered magnetic fields (e.g. in output files). We have verified 
that in this particular test problem the error for the BS versions (not reported here) are approximately 
10% less than for their UCT counterparts. 

The results for all UCT schemes are also plotted in Fig.O where the convergence rates are apparent. 
Note that in this kind of problems CENO-HLL is obviously by far the most performing code: the accuracy 
reached with a resolution of 32^ is comparable to that obtained by the second order methods with 128^. 

Concerning divergence-free properties, the particular settings of the problems are such that both 
[V • B]„„m = and [V • b]„tim = to within machine accuracy for the second order UCT schemes, 
because the invariance direction rj is made to coincide with the cell diagonal, a condition required by the 
double periodicity (independently on the value of the angle a). For CENO-HLL-UCT, typical values for 
the maximum values of [V • B]„„m and [V • h]num are 10~® and 10^^^, respectively. 

4.2 Rotated shock-tube problems 

These test problems involve the propagation of discontinuities defined by usual 1-D shock-tubes on a 
2-D computational plane, and are relevant to some main aspect considered here on the divergence-free 
condition. However, specific divergence-free properties may be hidden if one (or both) of the following 
special conditions hold: 

1. the initial magnetic field is uniform; 

2. the propagation direction lies along cell sides or diagonals. 

In the first case it is clear that any representation of [V • B]„„„i will be exactly zero for initial fields. 
Then its subsequent time evolution will only give a measure of the ability of a numerical scheme to 
preserve the initial [V • B]„ti„i, even if it is non- vanishing, while a characterizing aspect of any CT-based 
scheme is precisely the possibility to deal with discontinuous divergence-free fields. 
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Concerning the second case, the problem is the same as in the previous test, though here involves 
discontinuities and will be described in details. For initial symmetric conditions where all variables are 
defined as w(x,y) = w{£,), where ^ is a coordinate making an angle a with respect to the x axis, and 
they do not depend on the transverse i] coordinate, it is important to check at later times not only the 
evolution properties along the ^ coordinate but also the conservation of the transverse invariance. In 
particular, the V • B = condition expressed in the rotated coordinates is given by 

V • B = d^B^ + dr/B^ = 0, (58) 

and the equivalent condition B^{x,y) — const can be recovered only if (and all other variables, of 
course) are ry independent. However, in numerical applications based on standard Cartesian grids this 
condition may be achieved only if the ^ and rj directions are aligned with the cell diagonals (or the cell 
sides). In fact, any discontinuity front making a different angle will be discretized with unequal jumps 
along X and y. Thus, even when the [V • b]„„m ~ condition holds, the B^ = const relation does so only 
in approximate sense, with 0(1) jumps where discontinuities occur. More strictly related to the errors 
in the B^ variable, which has to be necessarily calculated from the interpolated Bi cell-centered fields, 
is the [V • BJnum variable, which in fact shows the highest jumps precisely at discontinuities. We have 
verified that, when the angle a — 45° is chosen, both the conditions B^ = const and [V • B]„ot„ = hold 
within machine accuracy, as in the previous test. 

The numerical domain for the oblique shock tube tests may be reduced, as cleverly suggested by 
Toth, to just a narrow strip [0,1] x [0, 2/7V], discretized with a x 2 grid (so that dx = dy). Shifted 
boundary conditions in the 77 direction are applied and a — tan^^ 2 w 63.4°. Each 2-D run, performed 
with N = 256, is compared with the corresponding 1-D test on a 1024 grid, by using the data along the 
X axis (the first row) at a final time tmax cos a (where ^max refers to the 1-D test). 

The initial left (L) and right (R) states of the three shock tube problems considered here are reported 
in Table @, and the final times are t^ax = 0.08 for ST-1, i^ax = 0.2 for ST-2, and i^ax = 0.1 for ST-3. 
ST-1 is a coplanar 2-D problem of converging shocks in an initially uniform magnetized background, ST- 
2 is a non-coplanar case involving Alfvenic discontinuities, and ST-3 is the famous (coplanar) Riemann 
problem involving the so-called compound (or intermediate) shock. Note that in ST-2 and ST-3 the 
magnetic field has jumps in the initial data and a ^ 45°, so none of the special cases indicated above 
apply. Other Riemann problems have been checked and the UCT schemes appear to behave well in all 
cases, including for example non-coplanar tests with compound shocks. The same three problems have 
been tested in LD for the CENO-LLF-UCT code on a symmetric N x N grid with a = 45°, while Toth 
shows just the first two tests; ST-1 with precisely the same settings and ST-2, the non-coplanar 2.5-D 
problem, with a — 45°. 

Here the most steepening slope limiters and a minimum amount of viscosity are used in our schemes: 
thus POSITIVE uses the Superbee (SB) limiter for all entries of L, cllf = 0.01, and the limiter in 
CENO-HLL is MC. In Table Q we report, for all tests and numerical schemes, the 5 average LI norm 
over the variables involved in each test, the LI norm of the errors in B^, the LI norms of variables 
[V • B]„„m and [V • b]„„m (respectively [V • B]avr and [V • b]avr), and their maximum absolute value 
over the computational domain (respectively [V • BJmax and [V • bjmax). The comparison between the x 
projection of the evolved quantities and the reference 1-D runs for our three UCT schemes are plotted in 
Fig. g. Fig. ®, and Fig. ®, for the ST-1, ST-2, and ST-3 tests, respectively. 

For problems involving sharp discontinuities, like shock tubes, schemes based on characteristics are 
clearly preferable since sharp limiters, which lead to more accurate results, may be often used there 
without producing spurious oscillations. This can be appreciated from the plots and may be measured 
more quantitatively from the reported errors in Table 0, which are the lowest for our POSITIVE scheme. 
If limiters less sharp than SB are used in POSITIVE-UCT the errors obviously increase: with MC we 
find S = 0.0146 for ST-1, S = 0.0153 for ST-2, and S = 0.0203 for ST-3; with the smoothest MM we find 
S = 0.0201 for ST-1, S = 0.0225 for ST-2, and S = 0.0298 for ST-3. Concerning the central schemes, the 
use of the MC limiter in the reconstruction step allows to produce accurate results with a rather low level 
of oscillations even in the absence of characteristics decomposition. These schemes are just less accurate 
at contact and Alfvenic discontinuities, since the related characteristic waves do not enter the HLL fiux 
definition. 

Finally, to appreciate the level of damage that monopoles can produce, we plot in Fig. the results 
of POSITIVE-BS for ST-3, where monopoles are free to arise (even in initial data since here none of the 
two special conditions apply) and grow. By also looking at the tabulated errors, it is clear that BS results 
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are systematically worst with respect to the correspondent UCT versions (see in particular the errors on 



4.3 The Orszag-Tang vortex problem 

A well-known model problem to study the transition to MHD turbulence is provided by the so-called 
Orszag-Tang vortex, which has been later adopted as a standard 2-D test for MHD shock-capturing codes. 
The initial conditions arc here Vx = — siny, Vy = sin a;, Bx = — siny, By = sin 2a;, p — 7^, p — j (so the 
the sound speed and the initial Mach number are both 1), Vz — Bz — 0. The computational domain is 
a square 0<x,y<2'KNxN box with periodic boundary conditions along both directions, while the 
output time is tmax — tt (note that in LD the magnetic field was normalized against Bq — l/\/47r). 

All our UCT schemes have been tested with the MC Hmiter (cllf = 0.1 in POSITIVE), at the 
resolutions of 50^, 100^, 200^, and 400^. In Table Q we report the averaged errors (over the six 
evolving variables). For each scheme, errors are measured with respect to a highest accuracy run at 
400^, obtained as the average of the results from the three UCT schemes. More qualitative comparisons 
can be appreciated in Fig. ((HJ, where gray-scale images of the temperature for the 200^ runs, for both the 
BS and UCT versions, are compared to the respective reference solution (the UCT run with 400^ grid 
points). Due to the symmetry of the problem, only the [0, tt] x [0, 2tt] first half of the domain is displayed. 
The first thing to be noticed is that the BS versions clearly produce incorrect results: the darkest, 
rather homogeneous features of the reference solutions appear more structured, with whiter filaments, 
and sometimes spurious oscillations are visible. By also reading Table it is then quite apparent 
that the three schemes behave similarly in this context. Due to the numerous steepened structures the 
order of accuracy is close to one in all cases, so that the errors are very similar. CENO-HLL-UCT is 
the most performing scheme, giving results almost identical to the Roe-type POSITIVE-UCT, while 
MC-HLL-UCT is of course more dissipative, but its results are still satisfactory. 

4.4 The fast rotor problem 

In a model problem to study the onset and propagation of strong torsional Alfven waves was presented 
and analyzed. A disk of radius r = 0.1 made up of dense fluid (rho = 10) rotates with high angular 
velocity {cj = 20) in a static, magnetized {B^ = b/y/in) background with uniform density and pressure 
{p = p — 1). The adiabatic index is 7 = 1.4. To conform with the final time \st = 0.15 and the same 
initial taper function is used (note that in LD the same problem was solved by the CENO-LLF-UCT code 
with t — 0.18 and no tapering). 

The numerical settings are identical to those employed in the previous test, and also the errors 
displayed in Table © are calculated in the same way. In Fig. © the magnetic pressure Pm = S^/2 is 
shown as isocontours diagrams for all schemes, first the BS and UCT runs at 200^ grid points, compared 
with the corresponding UCT reference run at 400^, which is also used to calculate the errors for the 
lower accuracy tests. Here the same remarks made above still apply: the accuracy order is low due 
to discontinuities, the two second order schemes behave very similarly (thus MC-HLL-UCT has to be 
preferred for its efficiency), while CENO-HLL-UCT gives the sharpest profiles, thus in this kind of model 
problems the accuracy in the reconstruction seems to be more important than the accuracy in resolving 
the Riemann structures. To conclude, we also notice that the BS schemes appear to behave correctly 
far from the rotating disk, where the waves are propagating outwards, while a lot of numerical noise is 
clearly formed inside the disk, where the numerical monopoles have time to accumulate, probably in a 
way similar to the inclined shock tube problems discussed in Sect. 4.2. 

5 Conclusions 

We have presented a method, first outlined in ^|, to construct Godunov-type schemes for the MHD sys- 
tem, named Upwind Constrained Transport (UCT). The main intent of our work is to assure that specific 
properties of the magnetic field, related to the basic divergence-free relation, enter as a built-in proper- 
ties also in the approximated systems. To that purpose, by taking advantage of the CT discretization 
technique, we have presented specific procedures to define consistent derivative approximations, recon- 
struction steps and approximate Riemann solvers all based on the staggered (or face centered) magnetic 
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field components hi chosen as primary data. A main advantage on this approach is that no cleaning 
procedures or ad hoc modifications of the form of the MHD conservative system are required. 
The main computational steps entering a UCT-based schemes are: 

1. reconstruction procedures based on the smoothness properties of the divergence-free B vector field, 
as represented in finite volume CT discretization (Sect. 2.2); 

2. the application of standard approximate Riemann solvers for the momentum and energy equa- 
tions, with the prescription that only variables not related to the divergence-free condition are 
reconstructed and participate to the upwind differentiation. As a benefit, among others, exact 
cancellation of numerical monopoles is assured (Sect. 2.3 and 2.4); 

3. a specific formulation of the approximate Riemann solvers for the induction equation (Sect. 2.3); 

4. a time integration procedure where no time-splitting is adopted. 

To demonstrate the validity and flexibility of our UCT method, we have finally applied it to a flux- 
limited Roe-type scheme (the positive scheme by Liu and Lax (14)), which proves to be accurate, robust 
and well-suited for more demanding applications requiring AMR techniques. This novel scheme has been 
then tested numerically on a standard set of model problems and compared to central-type second and 
third order schemes based on the two-speed HLL solver. 

We conclude by remarking that our method, defined here for the classical MHD system in regular 
structured grids, applies unchanged to the equations of special and general relativistic MHD (see TS'), and 
many procedures here presented may have a natural generalization for grid refinements and unstructured 
grids. 

Acknowledgments The authors thank G. Toth and another reviewer for their competent help in improving 
the manuscript. 



A Point-value formalism and third order procedures in the CENO- 



As shown in Sect. 2.3, to get high r > 3 order schemes in a finite-volume setting, besides a proper 
1 — th order reconstruction of W variables, a final averaging of the fi flux is needed. In the 3-D case the 
latter procedure is not cost-effective and more efficient implementations have been proposed for ENO- 
type schemes by Shu and Osher In this approach, point values u, instead of cell averages u, are 

advanced in time and flux values at the interfaces are directly reconstructed at the desired order by using 
flux point values fi(u). However, when applied to the MHD system (see |3|), this implementation is 
not suitable to take properly into account the divergence-free condition, since the flux reconstruction yet 
involves the Bi cell-centered values. Therefore, as in our previous work (LD), we use more appropriate 
flux reconstruction techniques. Let then consider the MHD equations, now in the form that comes out 
by applying the inverse operation of volume averaging to Eqs. and 



Here u = Uj^k,mi now cell-centered point values, constitute the new set of fluid primary data, and 
correspondingly the conservative flux two-point differences -Di(fi) are high order approximations of point- 
value first derivatives. In the induction equation, on the other hand, bi do not coincide with point-value 
representations of the staggered magnetic field components, those named bi which are needed in flux 
computations, since their volume averages must now return the sur/ace-averaged bi values. Therefore, bi 
components arc actually defined, now as primary data, in the same way fluxes are defined in Eq. H59|) . 
thus two-point differences Di{bi) give high order representations of point-value parallel first derivatives. 
Similarly, Dj{Ek) will be here high order representations of the staggered electric fields first derivatives, 
and their volume average must give back line-averaged electric fields. 



HLL-UCT scheme 




(59) 




(60) 
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It is important to notice that if hi components are evolved as primary data from Eq. (|65|1 . the point- 
value version of the divergence- free relation (|21|l . which is written as 

[V • hUm = D,ik) + Dy{by) + D,{b,) = 0, (61) 

will be preserved in time (if valid at i = 0) to within machine accuracy, exactly as in the finite-volume 
framework of Sect. 2.1. 

Our CENO-HLL-UCT scheme is then based on the following steps: 

1. For given cell-centered point values we reconstruct the face left-right point values using the 
appropriate CENO algorithm. The upwind fluxes fj^(wi,6i) and (w) are then evaluated by 
using HLL as in Sect. 2.3. 

2. For given values of these point- value representation of fluxes, the fj data are defined, for i — x,y, z 
and at the same points, as 

f^^{^~j^Di'\i,) (62) 

~ (2) 

where is a non-oscillatory approximation of the second derivative in the indicated coordinate. 
In this way, the difference Di{fi) provides a high order r > 3 accurate approximation of the flux 
first derivative. Concerning the magnetic fluxes, given the point values Ek at edge centers, the 
corresponding Ek data must be defined as 

4 = i?. - ^iWKE.) + Di'HE.)], (63) 
and similarly for x and y components. 

3. Eqs. (|59|1 and H6U|I can now be evolved in time by applying the Runge-Kutta algorithm of the 
appropriate (third) accuracy order. 

4. A final computational step is then needed to provide the relation between bi primary data and 
bi staggered point-value fields, those used in flux calculations. One has then to solve the implicit 
relations 

[I'^D^'\b,) = U, (64) 

for i — x,y, z, typically by means of iterative methods (see LD). As discussed in details in Sect. 2.4, 
the crucial point concerning how to measure magnetic monopoles is that the same algorithm to 
compute first derivatives for fiuxes should be actually applied to define the bi first derivatives 
in the [V • b]„„„i sum of Eq. H61() . which is exactly preserved only for hi primary data. Thus, 
when derivatives of [V • b]„„m are actually calculated starting from bi fields (obtained implicitly 
from Eq. 1)64(1 ). which is the definition actually relevant for numerical monopoles and that is in 
fact measured in our tests of Sect. 4, in a similar way as shown in Eq. (|62|l for fiuid fluxes, the 
[V • h]num variable will not be zero to within machine accuracy, though it can be made arbitrarily 
small depending on the accuracy of the inversion method employed. 

In the actual implementation of the code, if one wants to keep track of magnetic field-lines a slightly 
different approach, perfectly equivalent to that outlined above, can be followed. The induction equation 
(|^ may be substituted by 

^^[Ak{t)] = Ek, (65) 

where Ak is the point-value representation of the vector potential k component and Ek of the correspond- 
ing electric field component (see the line-averaged counterpart, Eq. IjlSfl ). The double non-oscillatory 
derivation in Eq. (|63|l must be now applied to, say, A^ rather than to E^'. 

A.=A,-^[d(^^{A^) + D<^HA^)], (66) 

and the divergence- free staggered 6j fields are now deflned (as in Eq. H14|l . proper for the finite- volume 
framework) as 

=^e«j,fci^j(ife), (67) 
to which the above remarks on [V • b]„t„„ apply unchanged. 
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Figure 1: The staggering of magnetic and electric vector fields in the CT framework. Only the 5+, 5*+ 
and cell faces arc visible. The arrows indicate the respective face normals, placed at intercell centers 
where 6j magnetic field components are defined, and the relative oriented contours for the application of 
Stokes' theorem, with arrows placed at edge centers where electric field components are defined. 
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Figure 2: The notation used for the reconstructed values inside a cell C'j^k,m (here a cut through the center, 
normal to the z direction, is shown) and the position of upwind fluxes, to be constructed via contribu- 
tions from neighbouring cells. Fluxes are either defined as two-state functions located at cell interfaces, 
{F^)j+i/2,k,m and {F'^)j,k+i/2,m, or as four-state functions located at cell edges, {'E^)j+i/2,k+i/2,m- 
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5ie S32 ^64 Si28 



POSITIVE-UCT 0.53097 0.12792 0.04273 0.01254 0.00322 
MC-HLL-UCT 0.60488 0.13133 0.04507 0.01392 0.00393 
CENO-HLL-UCT 0.31342 0.03759 0.00461 0.00056 0.00012 



Table 1: Averaged errors on the transverse velocity and magnetic field components for the oblique 2-D 

CP Alfvcn wave problem. The errors arc measured from the numerical solution at time t = I (one wave 
period), compared with the corresponding initial setting, for increasing resolutions and for the various 
UCT schemes. 




Figure 3: Convergence test for the oblique 2-D CP Alfven wave problem. Average Li errors on transverse 
T] and z wave components are shown in logarithmic scale for our three UCT schemes. Precise second 
order accuracy is achieved only asymptotically for POSITIVE-UCT and MC-HLL-UCT, due to the use 
of the MC limiter which tends to somehow sharpen wave profiles, while third order convergence is clearly 
reached by CENO-HLL-UCT, here employing MM in the reconstruction algorithm, already at very low 
resolutions. 
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1 


5/\/4^ 







Test ST-2: 


L 


1.08 


1.2 


0.01 


0.5 


0.95 




3.G/v^ 




Test ST-2: 


R 


1 











1 






2/V4^ 


Test ST-3: 


L 


1 











1 


0.75 


1 





Test ST-3: 


R 


0.125 











0.1 


0.75 


-1 






Table 2: Constant left (L) and right (R) states for the three oblique shock tube problems. 



Test ST-1 


(■) 




[V ■ B],,,,, 


[V ■ B],,-, 


[V ■ bjniux 


■ bjnvi 




POSITIVE-BS 


0.0225 


0.0051 


0.11 X 10^ 


0.30 X 10^ 








POSITIVE-UCT 


0.0126 


0.0019 


0.94 X 10^ 


0.26 X 10^ 


0.11 X 10-12 


0.67 X 10" 


15 


MC-HLL-BS 


0.0304 


0.0090 


0.12 X 10^ 


0.32 X 10^ 








MC-HLL-UCT 


0.0227 


0.0021 


0.10 X 10^ 


0.25 X 10^ 


0.11 X 10-^2 


0.44 X 10" 


15 


CENO-HLL-BS 


0.0320 


0.0092 


0.12 X 10-^ 


0.33 X 10^ 








CEXO-IILL-UCT 


0.0227 


0.0021 


U.ll X lO'"' 


0.27 X 


0.22 X 10"'> 


O.o7 X 10" 


-A 



Test ST-2 5 5B^ [V • B^x [V • B]avr [V • b]inax [V-b]avr 

POSITIVE-BS 0.0155 0.0016 0.65 x 10^ 0.34 x 10*^ '- - 

POSITIVE-UCT 0.0119 0.0006 0.70 x 10^ 0.29x10" 0.57 x lO'^^ 0.33 x 10"^^ 

MC-HLL-BS 0.0237 0.0025 0.61 x 10^ 0.35 x 10° 

MC-HLL-UCT 0.0200 0.0005 0.66 x 10^ 0.21 x 10° 0.57 x 10"^^ 0.22 x lO'^^ 

CENO-HLL-BS 0.0201 0.0025 0.71 x 10^ 0.41 x 10° 

CENO-HLL-UCT 0.0154 0.0006 0.63 x 10^ 0.25 x 10° 0.35 x 10"^ 0.14 x 10"^ 



Test ST-3 


5 


SB^ 




[V • B]avr 


[V • b]max 


[V • b]avr 


POSITIVE-BS 


0.0365 


0.0043 


0.16 X 10^ 


0.66 X 10° 






POSITIVE-UCT 


0.0165 


0.0005 


0.86 X 10^ 


0.30 X 10° 


0.14 X 10"" 


0.17 X 10-^^ 


MC-HLL-BS 


0.0586 


0.0064 


0.22 X 10^ 


0.13 X 10^ 






MC-HLL-UCT 


0.0295 


0.0003 


0.46 X 10^ 


0.19 X 10° 


0.28 X 10-13 


0.22 X 10-15 


CENO-HLL-BS 


0.0612 


0.0079 


0.25 X 10^ 


0.17 X 10^ 






CENO-HLL-UCT 


0.0257 


0.0004 


0.65 X 10^ 


0.23 X 10° 


0.13 X 10-6 


0.18 X 10-^ 



Table 3: Numerical errors for the three oblique shock tube problems. For the various schemes, the errors 
are calculated from the 2-D 256 x 2 run with respect to the corresponding 1-D high resolution run with 
1024 grid points. The displayed errors are the average 6 LI norm of all the involved variables, the error on 
the parallel field component (that supposed to remain constant), the maximum and averaged values 
of [V • B] 

ntini) and the maximum and averaged values of [V • b]„„m) which is not defined, of course, for 
non-CT schemes. 
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Figure 4: The oblique ST-1 shock tube problem. The numerical results from our three UCT schemes 
obtained on a 256 x 2 grid (symbols) are compared with the 1-D solution on a 1024 grid (solid line). 
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Figure 5: The oblique ST-2 shock tube problem. The numerical results from our three UCT schemes 
obtained on a 256 x 2 grid (symbols) are compared with the 1-D solution on a 1024 grid (solid line). 
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Figure 6: The oblique ST-3 shock tube problem. The numerical results from our three UCT schemes 
obtained on a 256 x 2 grid (symbols) are compared with the 1-D solution on a 1024 grid (solid line). 
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Figure 7: The oblique ST-3 shock tube problem, this time for the base scheme (BS) version of POSITIVE. 
Note the presence of errors induced by the numerical monopoles, here free to arise and grow in time. 



^50 Sion ^200_ 

POSITIVE-UCT 0.1601 0.0780 0.0296 
MC-HLL-UCT 0.1898 0.0920 0.0358 
CENO-HLL-UCT 0.1477 0.0657 0.0205 



Table 4: Averaged LI norms on the involved variables for the OT vortex problem at t = w. Errors are 
measured against a high resolution (400^) reference run. 



POSITIVE-UCT 0.1751 0.0888 0.0375 
MC-HLL-UCT 0.1770 0.0876 0.0331 
CENO-HLL-UCT 0.1525 0.0719 0.0239 



Table 5: Averaged LI norms on the involved variables for the fast rotor problem at f = 0.15. Errors are 
measured against a high resolution (400^) reference run. 
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Figure 8: Gray-scale images of the temperature T = p/ p distribution in the Orszag-Tang vortex problem. 
For each scheme, low resolution (200^) results for the BS and UCT versions are compared with the 
corresponding high resolution (400^) UCT reference run. 
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Figure 9: Contours of the magnetic pressure Pm = distribution in the fast rotor problem. For each 
scheme, low resolution (200^) results for the BS and UCT versions are compared with the corresponding 
high resolution (400^) UCT reference run. 
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